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ABSTRACT 


A  new  method  is  presented  for  the  analysis  of  MHD 
flow  in  the  entrance  region  of  a  parallel  plate  channel 
including  the  Hall  effect.  The  method  involves  a  linear¬ 
ization  of  the  inertia  terms  and  the  introduction  of  a 
stretched  axial  coordinate.  A  solution  is  obtained  for 
the  velocity  distribution  in  an  incompressible  fluid  which 
is  valid  from  the  duct  inlet  to  the  fully  developed  region. 

In  practice,  flows  in  which  the  Hall  effect  is  significant 
will  be  compressible.  Since  the  primary  objective  of  the 
present  analysis  is  to  develop  the  analytical  technique 
for  extending  the  linearization  scheme  to  MHD  flows  with 
Hall  effect,  the  constant  density  assumption  is  justifiable 
as  a  first  approximation.  From  this  solution  numerical 
results  are  obtained  for  the  entrance  length,  pressure 
distribution,  boundary  layer  thickness,  and  wall  shear  stress 
distribution  in  the  entrance  region  for  various  values  of 
the  Hartmann  number  and  the  Hall  parameter. 

Two  different  methods  of  analysis  of  fully  developed 
MHD  channel  flow  with  variable  conductivity  are  presented. 
Calculated  velocity  profiles  are  presented  by  the  two  methods 
for  several  conductivity  variations  across  the  channel. 

A  new  experimental  MHD  channel  flow  facility  designed 
to  use  potassium  chloride  (KC1)  as  the  working  fluid  is 
described.  Experimental  data  for  the  inlet  region  pressure 
distribution  in  this  facility  without  a  magnetic  field  are 
presented,  and  a  comparison  between  the  experimental  data 
and  analytical  prediction  is  made.  The  agreement  between 
the  experimentally  measured  entrance  region  pressure  drop 
and  the  theoretical  curve  is  considered  as  satisfactory 
verification  of  the  theoretical  model  employed  in  the 
present  investigation. 
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I 

INTRODUCTION 

In  recent  years,  two  promising  areas  of  engineering 
applications  of  magnetohydrodynamics  (MHD)  have  been  defined, 
namely  MHD  power  generation  and  MHD  acceleration.  A  recent 
survey  of  the  MHD  channel  flow  literature  by  Pai  (Ref.  1) 
indicates  that  most  of  the  analytical  studies  performed  to 
date  are  for  the  fully  developed  region,  the  region  far 
downstream  from  the  entrance  to  the  channel.  For  an  in¬ 
compressible  fully  developed  flow  at  constant  area,  flow 
conditions  are  independent  of  the  spatial  coordinate  in  the 
flow  direction. 

In  practice,  it  is  likely  that  a  significant  portion 
of  the  flow  in  MHD  generator  and  accelerator  channels  will 
occur  in  the  entrance  region  in  which  the  flow  undergoes  an 
adjustment  from  a  set  of  initial  conditions  to  the  fully 
developed  state.  In  addition,  the  practical  range  of  oper¬ 
ation  of  MHD  devices  indicates  that  the  Hall  coefficient  will 
be  significant  (Ref.  2). 

In  the  analytical  portion  of  the  present  investigation 
the  MHD  flow  in  the  entrance  region  of  a  parallel  plate 
channel  with  Hall  effect  present  is  analyzed.  A  linearization 
technique  developed  by  Sparrow,  Lin  and  Lundgren  (Ref.  3) 
for  non  MHD  flows  is  employed.  The  method  was  recently  ex¬ 
tended  to  MHD  flows  without  the  Hall  effect  by  Snyder  (Ref.  4). 

In  order  to  determine  the  validity  of  the  analytical 
solution,  experimental  data  are  presented  giving  measured 
velocity  profiles  and  pressure  distributions  in  the  entrance 
region  of  a  parallel  plate  channel.  The  working  fluid  for 
the  experiments  was  water.  A  comparison  between  analytical 
and  experimental  data  is  presented  which  verifies  the  validity 
of  the  analytical  method. 
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II 

PREVIOUS  WORK 


Relatively  few  studies  of  MHD  entrance  problems  are 
available  in  the  literature.  One  of  the  first  efforts 
made  to  estimate  the  MHD  entry  length  was  by  Sherd  iff 
(Ref.  5).  He  considered  a  linear  transient  problem  in 
which  a  flow  is  initiated  from  rest  in  an  infinitely  long 
channel  by  imposing  an  initial  condition  of  constant  mass 
flow  rate.  There  is  a  transient  period  in  which  the  flow 
settles  to  a  steady  state  flow  with  a  constant  pressure 
gradient.  The  time  required  for  the  flow  to  settle  to  a 
steady  state  condition  is  converted  to  an  entry  length  by 
multiplying  by  the  mean  velocity.  This  method  is  not 
precise;  as  acknowledged  by  Shercliff  himself,  and  even 
though  the  entry  length  may  be  given  to  the  correct  order 
of  magnitude  by  this  method,  the  details  of  the  flow  for 
the  transient  problem  would  not  be  expected  to  agree  with 
the  flow  in  the  actual  steady  state  entrance  problem. 

Roidt  and  Cess  (Ref.  6)  analyzed  the  entrance  problem  by 
applying  the  boundary  layer  technique  of  Schlichting  (Ref.  7). 
Since  nonsimilar  velocity  profiles  are  obtained  in  the 
boundary  layer  portion  of  the  flow,  a  series  solution  for 
the  stream  function  is  necessary.  The  labor  involved  in 
obtaining  higher  order  terms  in  the  series  becomes  excessive. 
Shohet,  Osterle,  and  Young  (Ref.  8)  obtained  a  numerical 
solution  to  the  entrance  problem  by  writing  the  governing 
equations  in  finite  difference  form  and  solving  them  on 
a  digital  computer.  The  same  problem  was  solved  also  by 
Hwang  and  Fan  (Ref.  9)  by  direct  numerical  integration  of 
the  governing  equations  in  finite  difference  form. 

All  of  the  above  analyses  used  the  simple  Ohm's  law 
in  the  form  J  =  a(E  +  v  x  B),  thus  excluding  the  Hall 
effect.  The  analysis  of  Chekmarev  (Ref.  10)  appears  to 
be  the  only  one  in  the  open  literature  in  which  an  effort 
is  made  to  include  the  Hall  current.  Chekmarev  considered 
an  approximate  linearized  version  of  the  entrance  problem 
which  will  be  subsequently  shown  to  be  a  special  case  of 
the  present  analysis.  The  principle  interest  of  Chekmarev 
was  in  the  formulation  of  the  problem  with  no  numerical 
results  given. 

In  the  present  analysis,  the  postulation  of  a  boundary 
layer  model  in  the  entrance  region  is  not  necessary.  A 
solution  for  the  velocity  distribution  is  obtained  which 
is  continuous  over  the  cross  sectional  area  of  the  duct 
from  inlet  to  the  fully  developed  region.  From  the  velocity 
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distribution  the  physical  quantities  of  interest  such  as 
the  entrance  length,  pressure  distribution,  wall  shear 
stress  distribution,  and  electric  current  distribution 
may  be  determined.  Because  of  the  presence  of  the  Hall 
current,  two  components  of  the  equation  of  motion  must  be 
retained  to  account  for  the  Hall  induced  transverse  velocity. 
Solutions  for  both  the  longitudinal  and  transverse  velocity 
components  are  obtained. 
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III 

THEORETICAL  INVESTIGATION 


3.1  FORMULATION  OF  THE  PROBLEM  AND  SOLUTION 
FOR  ARBITRARY  INITIAL  PROFILE 

The  geometry  of  the  MHD  parallel  plate  channel  is 
shown  in  Fig.  la.  The  relative  orientations  of  the  elec¬ 
tromagnetic  quantities  in  the  (x-y)  plane  are  shown  in 
Fig.  lb  for  generator  and  accelerator  modes  of  operation, 

A  constant  magnetic  field  Bz  is  applied  in  the  z-direction 
and  the  channel  is  assumed  to  be  infinite  in  the  y-direction. 
The  electric  field  E'  relative_to  the  moving  fluid  is 
related  to  the  electric  field  E  in  the  stationary  system 
by  the  relation 


E'  =  E  +  V  x  B 


(1) 


The  generalized  Ohm’s  law  including  ion  slip  and 
Hall  effect  has  been  discussed  by  Cowling  (Ref.  11)  and 
may  be  written  as 


J  - 


1+b 


(JOT  ,-r 
-  ( J  X 

(l+b)B 


B)  +  - - - 

B2 (1+b) 


(B  •  J)B 


(2) 


The  ion  slip  term,  b,  becomes  significant  only  for  cut  >  10. 
If  no  component  of  J  is  parallel  to  B,  the  term  J  •  B 
vanishes.  Assuming  the  magnetic  Prandtl  number  to  be  small 
relative  to  unity  implies  that  the  induced  magnetic  field, 
Bx,  is  small  compared  to  the  externally  applied  field,  Bz. 
With  these  assumptions,  the  x  and  y  components  of  J  may 
be  obtained  from  Eq.  (2)  as 


J 


x 


a 

(1+ft2) 


(3) 


a 

(i+a2) 


+  E^] 


(4) 
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where 


1+b 


£2 


CDT 

1+b 


The  governing  equations  of  motion  in  the  x  and  y 
directions  may  be  written  as 


P 


(u 

c)x  Sz 


ap  +  „  +  J  B 

Sx  bz2  ^ 


(5) 


p  (u 


dv 

dx 


+  w  — ) 
dz 


dP  d2v 

—  +  |x  — - 

dy  dz2 


J  B 


x  z 


(6) 


The  assumptions  implied  by  Eq.  (5)  and  Eq.  (6)  are 
as  follows: 


(1) 

(2) 


d  ,  ,r  dP 

- —  =  0  except  for  — 

dy  dy 


d  d 

Terms  involving  -r—  are  negligible  compared  to 

dx  dz 


(3)  The  contribution  to  the  body  force  by  the 
induced  field  Bx  is  negligible 

(4)  The  fluid  density  and  viscosity  are  assumed 
constant . 


Assumption  (1)  is  a  consequence  of  the  infinite  channel 
dimension  in  the  y-direction.  Assumption  (2)  is  a  standard 
boundary  layer  assumption  and  is  valid  except  in  the  im¬ 
mediate  vicinity  of  the  channel  entrance.  Assumption  (3) 
implies  a  small  magnetic  Reynolds  number,  a  valid 
assumption  for  most  laboratory  conditions  of  MHD  flows. 
Assumption  (4)  is  made  for  simplicity  because  the  primary 
objective  of  the  present  analysis  is  to  evaluate  the 
effect  of  the  Hall  current  on  the  entrance  region  flow 
development.  A  constant  property  assumption  is  reasonable 
as  a  first  approximation  to  study  the  qualitative  features 
of  the  flow  although  the  quantitative  results  would  certainly 
be  influenced  by  variable  properties. 

It  is  convenient  to  combine  Eq.  (5)  and  Eq.  (6)  into 
a  single  equation  by  introducing  complex  notation.  Multiplying 
Eq.  (6)  by  i  ~~\f~  1,  adding  the  result  to  Eq.  (5),  and 
substituting  for  the  current  components  from  Eq .  (3)  and 
Eq.  (4)  gives 
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dVc  bVc  dP  dP 

p(u  -  +  w  — -)  =  -  <—  +  i  — )  +  p  - 


dx  Sz 


dx  dy 


dz‘ 


aBz 

+  - —  [E  m-i)  +  E  (1+itt)  -  B  V  (1+ifi  )1 

x  y  z  c 

l4^2 


(7) 


where  Vc  =  u  +  iv  is  the  complex  velocity. 

In  (Ref.  4),  the  MHD  entrance  problem  without  Hall 
effect  is  analyzed  by  extending  the  linearization  technique 
of  Sparrow,  Lin  and  Lundgren  (Ref.  3)  to  the  MHD  case. 

Eq.  (7)  has  the  same  form  as  the  equation  of  motion  in 
the  absence  of  Hall  effect  and  thus  the  same  linearization 
technique  will  be  used.  Following  the  procedure  of  (Ref.  4), 
Eq.  (7)  is  linearized  to  the  form 


-  bVc 

e(x)U  -  =  A(x )  +  v 

dx 


+ 


oB. 


p(l+fi2) 


[Ex«l-i)  +  Ey(l+if2)  -  BzVc(l+ift)] 


(8) 


where  U  is  the  mean  velocity  defined  as 

h  h 

U  =  -  J*  V  dz  =  —  /  udz  (9) 

ii  c  h  . 

0  0 

The  mean  value  of  the  complex  velocity  V  is  equal 
to  the  mean  value  of  the  x-component  of  velocity  u  because 
of  the  physical  restraint  of  no  net  mass  flow  in  the  y- 

h 

direction,  i.e.  /  vdz  =  0.  Thus  U  is  a  real  quantity.  As 

0 

in  (Ref.  4),  e(x)  is  a  weighting  function  and  A(x)  includes 
the  complex  pressure  gradient  dP/dx  +  i  •  dP/dy  plus  the 
residual  of  the  inertia  terms  after  linearization.  In¬ 
tegrating  Eq.  (8)  over  the  upper  half  of  the  channel  cross 
section  determines  A(x).  Performing  the  integration  and 
utilizing  Eq.  (9)  gives 
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0  =  hA(x)  +  V 


3v, 


3z 


z=h 


crhB 


+ 


[E  (G-i)  +  E  ( 1+ifi )  -  B  U(l+ifl)] 


p(l+G2) 

Substituting  for  A(x)  from  Eq.  (10)  into  Eq.  (8)  gives 


(10) 


3v 

,,  3v 

eU  — —  =  • 

v. _ c 

3x 

h  3z 

z=h 


32V„  cB2  (1+ifl) 

+  v  - -  +  — - -  (U  -  V  ) 

3z2  p(l+ft2)  c 


(11) 


It  is  convenient  to  introduce  dimensionless  quantities 

V, 


a  =-  - 


z 

n  =  ~ 

U  U  U  b 


5  =  - 


0  =  -£. 


0B2h2 
2  -  z 


Ha  - 


M.(l+fi2) 


Hac2  =  (l+iG)Ha‘ 


R, 


=  £SL 


d3  = 


M- 


dx 

eR  h 
e 


(12) 


The  dimensionless  form  of  Eq.  (11)  becomes 


30  =  3f£  _  30 
3(3  3q2  3q 


q=l 


+  Hacz  (1  -  0) 


(13) 


with  boundary  conditions, 


39 


9 ( q~  1 )  =  0  (q=0)  =  0  0(0=0)  =  0  (q) 

3q 

Eq.  (13)  and.  associated  boundary  conditions  given 
by  Eq.  (14)  are  identical  to  Eq.  (9)  and  Eq.  (10)  of 
(Ref.  4)  and  thus  will  have  the  same  form  of  solution. 


(14) 
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The  solution  to  Eq .  (13)  satisfying  the  boundary  conditions 
of  Eq.  (14)  is 

q  _  Hac(cosh  Hacq  -  cosh  Hac) 
sinh  Hac  -  Hac  cosh  Hac 

Z-(72+Hac2)p 

Cne  (cos  7nq  -  cos  7r)  (15) 

n-1 

where  the  y  are  real  eigenvalues  and  are  the  roots  of  the 
transcendental  equation 


7n  =  tan  7n  n  =  1,  2,  . . .  (16) 

The  7n  are  the  same  eigenvalues  that  occurred  in  (Ref,  3) 
and  (Ref.  4),  and  the  first  twenty-five  eigenvalues  were 
tabulated  in  those  references.  During  the  course  of  the 
present  investigation,  a  computer  program  was  written  to 
evaluate  the  7n  since  more  than  twenty-five  terms  were  found 
to  be  necessary  to  obtain  convergence  under  certain  con¬ 
ditions.  After  these  calculations  were  performed,  a 
reference  by  Reisman  (Ref.  12)  was  found  in  which  the  first 
hundred  values  are  tabulated.  The  values  of  7n  are  tab¬ 
ulated  in  Table  1.  For  large  n,  the  7n  obey  the  recurrence 
relation, 


limit  (7  j  -  7n)  =  *  (17) 

n  ->  oo 


The  constants  Cn  are  complex  and  depend  on  the  initial 
velocity  profile  The  functions  fn(i})  =  cos7nq-cos7n 

form  an  orthogonal  set  over  the  interval  0  <  t)  <  1  which 
allows  determination  of  the  Cn  with  the  result 


(7^  +  Hac2)cos7n 

(Equation  continued) 
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1 

+  - - -  /  0O('])  (cosy  t)  -  cosy  )dq  0  8) 

sin2y 

'n  0 

The  first  term  of  Eq.  (15)  represents  the  fully 
developed  velocity  profile  in  complex  form.  It  is  seen 
to  be  identical  to  the  fully  developed  Hartmann  profile 
in  terms  of  the  complex  Hartmann  number  Hac .  That  the 
fully  developed  velocity  profile  with  Hall  effect  present 
could  be  put  in  the  form  of  the  Hartmann  profile  with  a 
complex  Hartmann  number  was  first  pointed  out  by  Yen 
(Ref.  13). 

It  should  be  emphasized  at  this  point  that  Eq.  (15) 
must  be  separated  into  real  and  imaginary  parts  to  be 
physically  meaningful.  Even  though  the  form  of  the  complex 
solution  is  identical  to  that  contained  in  (Ref.  4),  this 
apparent  simplicity  of  the  complex  solution  is  deceptive 
because  the  separation  of  Eq.  (15)  into  real  and  imaginary 
parts  is  a  rather  tedious  exercise  im  complex  algebra. 

This  separation  is  done  in  the  next  section  for  the  case 
of  uniform  initial  velocity  profile. 

3.2  VELOCITY  SOLUTION  FOR  UNIFORM  INITIAL  PROFILE 

For  a  uniform  initial  velocity  profile,  0o(q)  =  1  and 
Eq.  (18)  becomes 

C  =  - - -  (19) 

+  Hac2)cos7n 

The  real  part  of  Eq.  (15)_gives  the  dimensionless  longitudinal 
velocity  component  a  =  u/U_and  the  complex  part  gives  the 
transverse  velocity  6  =  v/U,  Because  of  the  lengthy  ex¬ 
pressions  involved,  it  is  convenient  to  list  the  fully 
developed  and  the  entrance  region  portions  of  the  velocity 
profile  separately.  The  first  term  of  Eq.  (15)  is  the 
fully  developed  profile  and  the  series  term  is  the  en¬ 
trance  region  profile.  Thus  the  dimensionless  velocity 
components  will  be  written 

a  =  —  =  a„  +  a  (20) 

U  1  e 


5  =  5  =  5 .  +  5p  (21) 

U 
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where  subscripts  f  and  e  refer  to  the  fully  developed  and 
entrance  portions  respectively.  After  a  considerable 
amount  of  complex  algebra  manipulations,  the  real  and 
imaginary  parts  of  Eq.  (15)  may  be  written  as 


cr,£  — 


coE„<t>  + 

C  2  +  D  2 
o  o 


(22) 


5f  = 


DoEo(^  ~  CoFo^> 

C  2  +  D  2 
o  o 


(23) 


CO  X  P 

nK 


a  =  2 


E 


n=l 


[cos7nT|  -  C0S7n]  [  X^^cos^  -  ?sin¥|3] 
[Xn2  +  ^2]cos  7n 


(24) 


5e  =  -  2 


00 

E 

n=l 


r^nP 


[cos7nq  -  coS7n]  [  Xnsinl,3  -  ’fcosE'(3] 


[X  2  +  ¥2]cos  7 


(25) 


n 


n 


C  -  sinh  A  cos  B  -  A  cosh  A  cos  B 

o  o  o  o  o  o 

+  B  sinh  A  sin  B  (26) 

o  o  o 

D  -  cosh  A  sin  B  -  A  sinh  A  sin  B 

o  o  o  o  o  o 

-  B  cosh  A  cos  B  (27) 

o  o  o 


Eo(rj)  =  Ao(cosh  AQr|  cos  BQq  -  cosh  Aq  cos  Bq) 

-  B  (sinh  A  n  sin  B  n 

o  o'  o  1 

-  sinh  Aq  sin  Bq)  (28) 
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FQ(r])  =  BQ(cosh  AQri  cos  Bqt]  -  cosh  Aq  cos  Bq) 

+  A  (sinh  A  ri  sin  B  n 
o  o  1  o  1 

-  sinh  A  sin  B  ) 
o  o 


r  =~\[l  +  SI2 


Hac  =  A  -  iB 
o  o 


n 


7/  +  Ha2 


¥  =  QHa2 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


The  longitudinal  and  transverse  velocity  components, 
given  by  Eq.  (22)  through  Eq,  (25)  may  be  used  to  calculate 
physical  quantities  of  interest  in  the  entrance  region. 
Before  these  calculations  are  presented  the  limiting  case 
of  cd  T  =0  will  be  considered  first.  From  the  definitions 
=  a)  t  /(I  +  b)  and  ¥  =  ftHa2,  it  is  clear  that  ¥  -  0  for 
cd  T  =0.  From  the  definitions  contained  in  Eq.  (26) 
through  Eq.  (35),  the  following  expressions  apply  for 
CD  T  =0. 


a  =  o 

(36) 

r  =  1 

(37) 
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B  =  0 
o 

A  =  Ha 
o 

C  =  sinh  Ha-  Ha  cosh  Ha 


Substituting 
through  Eq. 


+  2 


D  =  0 
o 


E  (q)  =  Ha(cosh  Harj  -  cosh  Ha) 


Vn>  “  0 


Hac  =  Ha 


Eq.  (36)  through  Eq.  (44)  into  Eq.  (22) 
(25)  gives 


_  Ha  (cosh  Hag  ~  cosh  Ha) 
sinh  Ha  -  Ha  cosh  Ha 


CO 

E 


n=l 


-(7n2  +  Ha2)p 


(cos  7nq  -  cos  7n) 


(7n2  +  Ha2 )cos  7n 


(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


(45) 


8  =  0 


(46) 


Eq.  (45)  is  identical  to  Eq.  (18)  of  (Ref.  4),  Eq.  (46) 
confirms  the  anticipated  result  that  the  transverse  velocity 
vanishes  for  <u  t  =  0. 

3.3  EVALUATION  OF  *(x) 

The  velocity  solution  given  by  Eq.  (22)  through  Eq.  (25) 
is  in  terms  of  the  stretched  axial  coordinate  0.  The 
physical  coordinate  x  is  related  to  the  stretched  coordinate 
0  from  Eq.  (12)  by  the  relation 
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P 


(47) 


Therefore  a  knowledge  of  the  variation  of  e  with  |3  is 
required  in  order  to  transform  from  the  (3  to  x  variable. 

In  (Ref.  4),  the  procedure  used  to  evaluate  e  was 
based  on  the  assumption  that  the  pressure  gradient  cal¬ 
culated  from  the  equation  of  motion  is  equal  to  the  pressure 
gradient  obtained  from  a  mechanical  energy  equation.  This 
procedure  will  be  used  in  the  present  analysis.  To  facil¬ 
itate  the  operations  which  are  to  follow,  it  is  convenient 
to  write  the  inertia  term  of  Eq.  (5)  in  the  form 


du  +  w  du  duu  duw 
dx  dz  dx  dz 


(48) 


Eq.  (48)  is  the  result  of  combining  the  continuity  equation 
with  the  inertia  term.  Combining  Eq.  (5)  with  Eq.  (4)  and 
Eq.  (48)  gives 

.  ,duu  .  duw^  dP  ,  ..  d2u 

dx  dz  dx  dz2 


CBZ 

+  - “T  [SHE  +  VB  )  +  (E  -  uB  )] 

1  +  ft2  X  Z  y 


(49) 


Integrating  Eq.  (49)  over  the  half  width  of  the  channel 
from  z  =  0  to  z  =  h  gives 


where 


u2dz  -  1  Su 
p  dx  h  dx  q  h  dz 


z=h 


aB  2U 

- - — ~  [^x  +  («v  ~  D] 

p(l  +  Q2)  x  y 


E 


X 


UB 


a>  = 


UB. 


(50) 

(51) 
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A  mechanical  energy  equation  may  be  obtained  by 
multiplying  Eq.  (5)  by  u  with  the  result 


1_  3uu2  1^  3wu2  __  _  u  3P  +  vu  32u 
2  3x  2  3z  p  3x  3z2 

aB 

+  - 5 — -  [Q(uEt  4-  uvB  )  +  (uE  -  u?B  )  ]  (52) 

p(l  +  fi2)  X  2  y  2 


Integrating  Eq .  (52)  gives 


1  3P 
p  3x 


_3_ 
hU  3x 


dz  + 


/  (SH)=dZ 

hU  Q  3z 


gB22^ 
pCi  +  a2) 


h 

/  uvdz) 
0 


($  -  -i-  /  U2dz)] 


hU‘ 


(53) 


Equating  the  pressure  gradient  expression  from  Eq.  (50) 
and  Eq.  (53)  and  introducing  dimensionless  quantities  gives 
the  following  expression  for  e. 


e(Ha,Q,p)  = 


2-  f  ( a2  -  *d)«lT, 

33  0  _ 2 


(54) 


5  Ti=l+Ha2^  (a2-l)dr\-Cl  /  o6drj] 


3a 


W 


0 


The  right  side  of  Eq.  (54)  is  a  known  function  of  3, 
Ha,  and  from  Eq.  (22)  through  Eq.  (25).  The  physical 
coordinate  x  is  then  given  by  Eq.  (47). 

Eq.  (54)  shows  that  €,  the  weighting  function,  is  a 
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function  of  (3  as  well  as  the  electromagnetic  parameters  Ha 
and  £1.  The  integrals  appearing  in  Eq .  (54)  were  evaluated 
numerically.  The  relation  between  the  physical  coordinate 
x  and  the  stretched  coordinate  p  was  determined  from  Eq . 

(47)  which  involves  a  numerical  integration  of  e  with  respect 
to  p.  The  results  are  shown  in  Figs.  2,  3,  and  4.  Shown 
also  for  comparison  are  the  curves  for  e  =  1,  an  assumption 
made  by  Chekmarev  (Ref.  10).  The  assumption  e  =  1  is  seen 
to  be  a  poor  one  as  Ha  increases  at  cot  =  0.  For  cot  =  10, 
the  e  =  1  curve  is  approached  more  closely.  Basing  the 
entrance  length  determination  on  the  assumption  e  =  1  results 
in  a  longer  entrance  length  prediction  than  that  based  on 
a  variable  e  given  by  Eq.  (54). 

3.4  THE  ENTRANCE  LENGTH 

For  incompressible,  constant  area  flow,  a  fully 
developed  velocity  profile  is  approached  which  is  independent 
of  the  axial  coordinate.  Because  the  fully  developed  pro¬ 
file  is  approached  asymptotically,  the  definition  of  the 
entrance  length  is  arbitrary.  The  customary  definition  of 
the  entrance  length  is  the  distance  required  for  the  center- 
line  velocity  to  become  a  specified  fraction  of  the  fully 
developed  centerline  velocity.  The  difference  between  the 
fully  developed  centerline  velocity  and  the  centerline 
velocity  at  any  position  p  may  be  written  from  Eq.  (24)  as 

oo  -X  p 

Ze  n  (cost  — 1 )  ( X  cosfp  -  fsin'ip) 

- II - - -  (55) 

(X  2  +  lf2)cos  7 
n  'n 

n=l 


The  entrance  length  will  be  defined  as  that  distance 
at  which  a_f(r]  =  0)  -  a(q  -  0)  =  0.05.  This  definition 
corresponds  to  the  location  at  which  the  centerline  velocity 
is  95%  of  the  fully  developed  centerline  value.  Entrance 
length  values  are  shown  in  Fig,  5  for  Hartmann  number 
values  of  Ha  =  2,  4,  6  and  in  Fig.  6  for  values  Ha  =8,  10. 
For  a  given  value  of  Ha,  the  entrance  length  is  seen  to 
increase  with  increasing  values  of  cdt.  For  the  smaller 
values  of  Ha,  the  entrance  length  values  appear  to  be 
approaching  an  asymptotic  limit  as  cot  increases.  The 
asymptotic  limit  is  not  so  obvious  for  the  larger  values 
Ha  -  8  and  Ha  -  10  since  cot  =  10  was  the  largest  value 
used  in  the  calculations.  However,  an  examination  of  the 
complete  equation  of  motion  given  by  Eq.  (7)  shows  that 
in  the  limit  as  cot  oo  with  B  and  E  held  constant,  the 
electromagnetic  body  force  term  vanishes.  Thus  as  cot  ->  oo, 
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the  entrance  region  flow  approaches  the  non  MHD  entrance 
region  as  a  limit.  Additional  justification  of  this  con¬ 
clusion  is  provided  by  an  examination  of  the  influence  of 
cox  on  the  fully  developed  velocity  profiles. 

3.5  INFLUENCE  OF  <yrON  FULLY  DEVELOPED  VELOCITY  PROFILE 

Fully  developed  velocity  profiles  are  shown  in  Figs. 
7,  8,  9  for  Ha  values  of  2,  6,  10  and  oyr  values  of  0  and 
10.  Shown  also  for  comparison  is  the  profile  for  fully 
developed  flow  in  the  absense  of  a  magnetic  field,  the 
well  known  parabolic  distribution  given  by  the  relation 


-  =  -  [1  -  (-)2]  (56) 

U  2  h 


As  out  increases,  the  curves  approach  the  parabolic  dis¬ 
tribution.  Thus  the  parameters  Ha  and  cot  distort  the 
velocity  profiles  in  opposite  trends.  Increasing  Ha  tends 
to  flatten  the  profile  in  the  center  of  the  channel,  re¬ 
sulting  in  an  increased  shear  stress  at  the  wall,  whereas 
increasing  art  causes  the  profile  to  become  more  nearly 
parabolic  in  shape  ,  resulting  in  a  decrease  in  the  wall 
shear  stress.  This  trend  in  the  fully  developed  region 
suggests  that  in  the  entrance  region,  the  boundary  layer 
thickness  should  increase  as  oyr  increases. 

3.6  INFLUENCE  OF  at  ON  ENTRANCE  REGION  BOUNDARY  LAYER  THICKNESS 


The  definition  of  the  boundary  layer  thickness  in 
either  external  or  internal  flows  is  arbitrary  because  of 
the  non-uniqueness  of  the  boundary  layer  thickness.  In 
internal  flows,  a  frequently  used  definition  is  the  distance 
from  the  wall  at  which  the  velocity  becomes  a  specified 
fraction  of  the  centerline  velocity.  In  the  present 
analysis,  the  boundary  layer  thickness  will  be  defined  as 
the  distance  from  the  wall  at  which  the  local  velocity  is 
90%  of  the  centerline  velocity.  This  definition,  although 
arbitrary,  will  serve  the  purpose  of  illustrating  the 
influence  of  oyr  on  the  boundary  layer  thickness. 

The  entrance  region  boundary  layer  thickness  is  shown 
plotted  in  Figs.  10,  11,  12  for  Ha  values  of  2,  6,  10  and 
cox  values  of  0  and  10.  The  curves  are  terminated  at  the 
end  of  the  entrance  region  that  was  previously  defined  as 
the  position  at  which  the  centerline  velocity  is  95%  of 
the  fully  developed  centerline  value.  The  boundary  layer 
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thickness  so  defined  is  seen  to  increase  as  oyr  increases, 
a  result  qualitatively  compatible  with  the  observed  effect 
of  cot  on  the  fully  developed  velocity  profiles. 

3.7  ENTRANCE  REGION  PRESSURE  DISTRIBUTION 

The  pressure  distribution  in  the  entrance  region  may 
be  obtained  by  integrating  either  Eq.  (50)  or  Eq .  (53) 
from  the  entrance,  x  =  0,  to  any  position  x.  Because  of 
the  simpler  form  of  Eq.  (50),  it  will  be  used  to  determine 
the  pressure  distribution.  To  emphasize  the  entrance  region 
effect  on  the  pressure  distribution,  it  is  convenient  to 
separate  the  pressure  drop  into  a  fully  developed  portion 
and  an  entrance  region  portion.  This  separation  can  be 
accomplished  by  using  the  identity 


du 

dz 


z=h 


duf 

dz 


+ 


z=h 


d(u  -  uf) 
dz 


z=h 


(57) 


where  Uf  is  the  fully  developed  velocity.  Substituting 
Eq.  (57)  into  (50)  and  integrating  from  x  =  0  to  x  =  x 
gives 


P  -  p  -|  .  ,,  r  du, 

-2 - ,  i  /  u2dz  -  U2  -  ^  /  — £ 

p  h  Q  11  0  Sz 


dx 


z-h 


x 


V  f  d(u  -  V 

h  „  dz 


a  B  2Ux 

dx - 2 -  [ft*  +  (*  -  D]  (58) 

z=h  p(l  +  ft2)  X  7 


Writing  Eq,  (58)  in  dimensionless  form  gives 


Po  -  p 


!pir 


—  (-) 
Re  h 


/  da. 


ri=l 


+  Ha[ft*x  +  (*  -  1)] 


}  P  d(a  -  a,) 

] 

+  2 

J  a2dq  -  1  -  J  e 

dp 

0  0 

rrl  | 

(59) 
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The  first  bracketed  term  of  Eq.  (59)  is  the  pressure 
drop  resulting  from  the  fully  developed  flow  and  the  second 
bracketed  term  is  an  entrance  region  correction  term.  Re¬ 
writing  Eq,  (59)  in  terms  of  a  fully  developed  friction 
factor  and  an  entrance  region  correction  term  gives 


P  -  P  , 

°  -  f  —  (-)  +  K 


ipU2 


Re  h 


(60) 


where 


f  =  _  2  |  — 
dt) 


rpl 


+  Ha2[ft$x  +  ($  -  1)] 


(61) 


K  =  2  (  /  -  1 

0 


3 


- / 


3(a  -  ou) 


0 


dri 


1=1 


dp  (62) 


The  fully  developed  friction  factor,  f,  is  a  function 
of  both  components  of  the  electric  field  3^  and  3>y.  The 
entrance  region  pressure  correction,  K,  is  independent  of 
the  electric  field,  however.  It  is  interesting  to  observe 
that  the  Hall  current  induced  transverse  velocity  component, 
5,  does  not  appear  in  the  longitudinal  pressure  distribution. 
It  would  appear  in  the  transverse  pressure  distribution, 
however,  since  the  transverse  velocity  induces  a  transverse 
pressure  distribution. 

The  entrance  region  pressure  correction  is  shown  in 
Figs.  13,  14,  15  for  Ha  values  of  2,  6,  10  and  cut  values 
of  0  and  10.  An  examination  of  the  curves  shows  that  for 
a  fixed  value  of  Ha,  K  increases  with  increasing  cor.  The 
influence  of  oyr  becomes  more  pronounced  as  Ha  increases. 

For  a  fixed  value  of  cot,  K  decreases  with  increasing  Ha. 

The  curves  are  terminated  at  the  end  of  the  entrance  region. 

For  the  fully  developed  region,  the  pressure  gradient 
becomes  constant  and  is  given  by  the  fully  developed 
friction  factor  defined  by  Eq.  (61).  Thus  it  is  clear 
that  the  value  of  K  should  asymptotically  approach  a  con¬ 
stant  value  in  the  limit  as  x  ->  «?.  This  conclusion  is 
also  evident  from  the  definition  of  K  given  by  Eq.  (62). 

The  asymptotic  values  of  K  are  shown  on  Figs.  13,  14,  15. 
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Since  the  curves  were  terminated  at  the  arbitrarily  defined 
entrance  length,  i,e,  the  position  at  which  the  centerline 
velocity  becomes  95%  of  the  fully  developed  centerline 
velocity,  the  asymptotic  value  is  not  reached  by  the  curves 
The  total  pressure  in  the  entrance  region  should  be  based 
on  the  asymptotic  value  of  K.  The  variation  of  the  asymp¬ 
totic  value  of  K  with  cut  is  shown  in  Fig.  16  for  Ha  values 
of  2,  6,  10. 

3.8  ENTRANCE  REGION  WALL  SHEAR  STRESS  DISTRIBUTION 

The  shear  stress  at  the  wall  is  given  by  the 
expression 


5u  _  joU  5a  [ 

5z  z=h  h  5t||t)=1 


(63) 


The  negative  sign  is  used  to  make  t  a  positive  quantity 
since  the  velocity  gradient  at  the  $all  is  negative.  De¬ 
fining  a  dimensionless  shear  stress  as 


T*  =  — S—  (64) 

<n&> 

h 


and  evaluating  the  velocity  gradient  from  Eq.  (20),  (22), 

and  (24)  gives 


T*  = 
W 


[ sinh  Ao  cos  Bo][Co(Boa  -  AQS)  -  2AqBoDo] 


C  2  +  D  2 
o  o 


+ 


f cosh  A  sin  B  ][D  (B2  -A2)  +2ABC  ] 
L  o  oJLo  o  o'  oooJ 


C  2  +  D  2 
o  o 


co  -x  p 

n  2 


+  2 


E 


n=l 


e  7  cos  Tp  -  Y  sin  ¥0] 


[Xn2  +  T2] 


(65) 
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The  first  two  terms  in  Eq.  (65)  give  the  wall  shear  stress 
in  the  fully  developed  region  and  the  summation  represents 
the  contribution  due  to  the  entrance  region. 

The  variation  of  the  wall  shear  stress  in  the  entrance 
region  is  shown  in  Figs.  17,  18,  and  19.  The  results  are 
plotted  as  a  ratio  of  local  wall  shear  stress  to  the  fully 
developed  value.  Since  the  initial  velocity  profile  at  the 
channel  entrance  was  assumed  to  be  uniform,  a  discontinuity 
in  the  velocity  distribution  at  the  wall  occurs,  resulting 
in  an  infinite  value  of  shear  stress  at  the  inlet.  This  is 
not  a  physically  meaningful  result  but  is  a  consequence  of 
the  assumption  of  an  initial  uniform  profile,  a  physically 
impossible  condition  for  a  viscous  fluid. 


3.9  INFLUENCE  OF  INITIAL  VELOCITY  PROFILE  ON  ENTRANCE  LENGTH 

The  solution  for  the  entrance  region  velocity  profile 
is  given  by  Eq.  (15)  with  the  constants  Cn  given  by  Eq. 

(18).  The  0o(q)  appearing  in  Eq.  (18)  is  the  initial 
velocity  profile  at  the  channel  entrance.  The  results 
presented  so  far  have  been  for  uniform  initial  profile, 
i.e.  0o  =  1.  In  an  actual  MHD  channel,  the  initial  pro¬ 
file  would  not  be  uniform. 

In  order  to  assess  the  influence  of  the  initial 
velocity  profile  on  flow  development  in  the  entrance  region, 
a  family  of  initial  profiles  was  assumed  described  by  the 
relation 


0OW  =  ao(l^  =  ~-+i  -  t1  -  tin]  (66> 


where  N  is  a  constant.  This  profile  is  seen  to  satisfy 
the  relation 


1 

J  S0(n)dT!  =  1  (67) 

0 


which  is  required  by  the  definition  of  mean  velocity  and 
the  condition  of  no  net  mass  flow  in  the  transverse  direction. 
This  family  of  initial  profiles  was  chosen  because  it  allows 
representation  of  the  laminar  fully  developed  non-MHD  flow 
as  well  as  flows  with  steeper  velocity  gradients  near  the 
wall  than  the  laminar  fully  developed  gradient.  Thus  a 
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turbulent  type  profile  could  be  approximated  by  taking  N 
sufficiently  large.  The  uniform  initial  profile  is  also 
included  as  a  limiting  case  in  the  limit  N  -*■ 

Numerical  results  were  obtained  for  the  three  values 
N  =  2,  6,  10  with  cox  =  0.  This  latter  condition  implies 
that  the  transverse  velocity  component  vanishes  and  thus 
9  -  a.  The  entrance  length  variation  with  N  and  Ha  is 
shown  in  Fig.  20.  For  N  =  2,  corresponding  to  a  parabolic 
initial  velocity  profile,  the  variation  of  entrance  length 
with  Ha  is  seen  to  begin  at  0  for  Ha  =  0,  reach  a  maximum 
value  at  Ha  =  4,  and  approach  asymptotically  a  value  in¬ 
dependent  of  Ha  for  Ha  >  10.  This  qualitative  feature  of 
a  maximum  value  of  entrance  length  occurring  is  present 
only  for  the  case  N  =  2  and  can  be  explained  as  follows. 

For  a  non-MHD  flow,  N  =  2  corresponds  to  the  fully  developed 
laminar  profile.  Thus  for  Ha  =  0,  N  =  2,  the  inlet  flow  is 
already  fully  developed,  giving  an  entrance  length  of  zero. 

As  N  increases,  the  entrance  length  variation  with  Ha  is 
seen  to  approach  the  curve  corresponding  to  a  uniform 
initial  profile. 

A  measure  of  the  influence  of  the  initial  velocity 
profile  on  the  entrance  region  flow  development  can  be 
obtained  by  comparing  the  wall  shear  stress  of  the  initial 
profile  with  the  fully  developed  wall  shear  stress  for  a 
given  value  of  Hartmann  number.  From  Eq.  (66)  the  dimension¬ 
less  shear  stress  at  the  wall  is  given  by 

C  -  N  +  1  t68> 


The  fully  developed  MHD  wall  shear  stress  for  cdt  =  0  is 
given  by 


T 


* 

W 


_ Ha2 

Ha  coth  Ha  -  1 


(69) 


For  large  Ha,  the  difference  between  Eq.  (69)  and  Eq.  (68) 
may  be  written 


N 


(70) 


Eq.  (70)  shows  that  for  large  values  of  Ha,  the  fully 
developed  wall  shear  stress  and  the  initial  wall  shear 
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stress  become  equal  when  Ha  =  N.  This  conclusion  is 
illustrated  in  Fig.  21  in  which  the  entrance  region  center- 
line  velocity  is  plotted  for  N  =  2,  6,  10  and  Ha  =  2,  6, 

10*  From  the  curves  it  is  clear  that  when  Ha  =  N,  there  is 
very  little  change  in  the  centerline  velocity.  Since  the 
primary  mechanism  involved  in  adjusting  an  initial  velocity 
profile  to  a  fully  developed  profile  is  the  diffusion  of 
vorticity  generated  at  the  wall  into  the  central  portion  of 
the  flow,  if  the  inlet  wall  shear  stress  is  equal  to  the 
fully  developed  wall  shear  stress,  there  will  be  no  signif¬ 
icant  adjustment  of  vorticity  necessary  and  hence  no  velocity 
entrance  length.  This  conclusion  would  not  be  valid  for 
a  compressible  flow  with  variable  properties,  however. 

3.10  A  GENERALIZED  MOMENTUM  INTEGRAL  ANALYSIS 

The  previous  analysis  of  the  MHD  entrance  region 
problem  did  not  use  a  boundary  layer  -  core  flow  model. 
Instead  the  equation  of  motion  was  linearized  in  such  a 
way  that  a  solution  for  the  velocity  profile  was  obtained 
which  is  continuous  from  the  entrance  of  the  channel  to 
the  fully  developed  region.  If  a  boundary  layer  -  core 
flow  model  is  used,  two  techniques  are  available,  (1)  the 
momentum  integral  method  and  (2)  a  non  similar  series 
expansion. 

Several  investigations  of  the  MHD  entrance  region 
problem  based  on  the  momentum  integral  method  have  appeared 
in  the  literature  in  recent  years,  Moffatt  (Ref.  14)  was 
the  first  investigator  to  apply  the  momentum  integral 
method  to  the  study  of  MHD  boundary  layer  formation  over 
a  flat  plate.  He  assumed  a  parabolic  velocity  distribution 
in  the  boundary  layer. 

Maciulaitis  and  Loeffler  (Ref,  15)  applied  the  momentum 
integral  method  and  got  a  closed  form  solution  by  assuming 
the  velocity  distribution  to  have  a  parabolic  variation 
and  the  centerline  velocity  to  vary  down  the  channel. 

Tan  (Ref.  16)  also  found  a  closed  form  solution  by 
using  the  momentum  integral  method  and  assuming  a  Hartmann- 
like  velocity  distribution  previously  used  by  Moffatt 
(Ref.  14)  who  had  solved  the  problem  numerically. 

Dhanak  (Ref.  17)  applied  the  classical  Karman-Pohlhausen 
method  to  the  entrance  region  and  solved  the  resulting 
equation  by  an  iterative  scheme. 

Although  several  authors  have  applied  momentum  integral 
methods  to  MHD  entrance  flow,  using  simple  parabolic, 
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Hartmann-like  or  Karman-Pohlhausen  velocity  distributions, 
these  works,  rather  than  establish  the  validity  of  the 
momentum  integral  method.,  have  created  controversy  as  to 
whether  or  not  one  can  validly  apply  such  a  method  to  MHD 
flows,  Hugelman  and  Haworth  (Ref.  18)  reasoned  that  the 
classical  Karman-Pohlhausen  method  does  not  take  into 
account  enough  boundary  conditions.  They  suggested  using 
a  fifth-degree  polynomial.  Heywood  and  Moffatt  (Ref.  19) 
state  that  integral  method  results  are  too  much  dependent 
upon  the  velocity  profile  assumed  and  that  Dhanak’s  solution, 
in  particular,  did  not  proceed  far  enough  downstream  to 
show  clearly  the  effect  of  the  magnetic  field. 

One  of  the  most  severe  limitations  of  the  classical 
momentum  integral  method  in  which  fourth-degree  polynomial 
profiles  are  assumed  is  that  the  solution  cannot  be 
obtained  over  the  entire  range  from  entrance  to  the  fully 
developed  region.  This  is  a  limitation  for  an  MIID  entrance 
region  flow  as  well. 

In  the  present  investigation  the  momentum  integral 
method  has  been  generalized  to  include  a  polynomial  velocity 
profile  of  arbitrary  degree.  It  is  found  that  increasing 
the  degree  of  the  assumed  polynomial  velocity  profile  per¬ 
mits  a  solution  which  is  continuous  from  the  entrance  region 
to  the  fully  developed  region. 

The  method  will  be  developed  for  the  case  cot  -  0 
which  means  that  only  the  x-component  of  the  equation 
of  motion  is  needed.  Eq.  (5)  with  cur  =  0  becomes 

p(u  —  +  w  =  -  —  +  n  Sfi*  +  OB  (E  -  llBJ  (71) 

dx  dz  dx  dz2  z  y  z 

with  the  continuity  equation  given  by 


du  [  dw 
dx  dz 


0 


(72) 


Eq.  (71)  is  integrated  over  the  boundary  layer  thickness 
A  to  obtain  the  momentum  integral  equation  in  the  form 


d£* 

dx 


2@*aB  2 
_ z _ 

PUs 


o*  c*  dU  aB 

2-  (2  +$-)(—£+  -2 

0*  dx  p 


(73) 
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where  5*  and  0*  are  the  displacement  thickness  and 
momentum  thickness  respectively,  defined  as 

A 

(74) 


(75) 


The  velocity  us  is  the  uniform  core  velocity  outside_the 
boundary  layer.  A  relationship  between  &*,  us,  and  U, 
the  average  velocity  may  be  obtained  from  a  mass  balance 
given  by  the  expression 


5*  -  f  (1  ~  — )dz 
0  us 

A 

e*  =  /  —  (i  -  -±-)dz 

n  U„ 


A  h 

Uh  =  /  udz  +  J  u  dz  (76) 

b 

0 

Adding  and  subtracting  the  term  f  ugdz  to  the  right  side  of 
Eq.  (76)  and  rearranging  gives  0 


(77) 


To  solve  Eq.  (73)  requires  an  assumed  velocity  profile. 
The  conventionally  used  profile,  following  the  classic  work 
of  von  Karman  and  Pohlhausen,  is  a  polynomial  of  degree  four. 
The  fourth-degree  polynomial  profile  has  certain  limitations, 
however,  the  most  serious  being  that  a  solution  cannot  be 
obtained  which  is  valid  from  the  entrance  region  to  the 
fully  developed  region.  At  a  certain  distance  downstream 
from  the  entrance,  a  velocity  overshoot  occurs  in  the 
boundary  layer,  that  is  the  local  boundary  layer  velocity 
becomes  larger  than  the  local  free  stream  velocity.  Although 
such  an  overshoot  is  possible  in  a  compressible,  variable 
property  fluid,  there  appears  to  be  no  physical  mechanism 
which  would  allow  a  velocity  overshoot  in  an  incompressible, 
constant  property  fluid.  Thus  the  momentum  integral  solution 
is  considered  to  be  valid  only  to  the  position  at  which 
boundary  layer  overshoot  occurs. 


24 


A  E  DC-T  R  -67  *79 


In  the  present  investigation,  the  momentum  integral 
method  is  generalized  to  allow  a  polynomial  velocity 
profile  of  arbitrary  degree.  It  is  shown  that  for  a 
given  Hartmann  number,  taking  a  polynomial  profile  of 
sufficiently  high  degree  will  always  allow  a  solution 
which  is  valid  over  the  entire  entrance  region. 


Consider  a  boundary  layer  profile  expressed  as  a 
polynomial  of  degree  N  of  the  form 


N 


u 


u 

s 


£ 

n=1 


n 


(78) 


where 


z 

A 


The  coefficients  an  will  be  functions  of  x  and  N  boundary 
conditions  are  required  to  determine  the  N  values  of  an. 
In  the  present  analysis,  three  boundary  conditions  are 
applied  at  the  wall,  -q*  =  0,  and  the  remaining  N  -  3  are 
applied  at  the  outer  edge  of  the  boundary  layer,  r)*  =  1. 
In  particular  the  boundary  conditions  become 


at  q*  =  0 


u  -  0  v  =  0 


v  ifn  =  1  dP 

dz2  p  dx 


CB  E 

z  y 


P 


dus 

u  ( - -  + 

s  dx 


(79) 


at  t]*  =  1 


u 


-  u 

s 


du  _  52u 
bz  bz2 


c>zm 


0 


m  =  N  -  3 


(80) 


The  first  two  conditions  of  Eq.  (79)  are  the  velocity 
non  slip  conditions  and  the  third  condition  is  the  equation 
of  motion  evaluated  at  the  wall.  The  first  condition  of 
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Eq.  (80)  requires  the  boundary  layer  velocity  to  match  the 
free  stream  velocity  at  the  outer  edge  of  the  boundary  layer 
and  the  remaining  conditions  require  the  first  N  -  3  de¬ 
rivatives  of  the  velocity  profile  to  vanish  at  the  outer 
edge . 


Substituting  the  boundary  conditions  from  Eqs.  (79) 
and  (80)  into  Eq.  (78)  gives 


a  -  -  -  (81) 

2  2 

N 

Y  an  =  1  (82) 

n=l 


N-(m-l) 


Z 

n— 1 


[  (n+m-1) (n+m-2) 


n  ]a 


n+m-1 


0  m-1,2, 


N-3  (83) 


where 


X* 


v  dx 


(84) 


The  parameter  X*  defined  by  Eq.  (84)  is  the  boundary 
layer  shape  factor  and  it  is  clear  that  each  of  the  a_ 
coefficients  will  be  a  function  of  X*  and  thus  a  function 
of  x. 


Eqs.  (81)  through  (83)  are  a  set  of  N  equations  for 
the  N  unknown  an  coefficients.  For  large  N,  the  direct 
solution  of  the  equations  in  explicit  form  is  not  feasible. 

By  starting  with  a  different  form  for  the  polynomial  velocity 
profile,  however,  it  is  possible  to  obtain  an  exact  re¬ 
currence  relation  for  all  the  an  in  terms  of  X*,  valid  for 
arbitrary  degree  N. 

For  this  purpose,  consider  the  function 

—  =  (1  -  ti*)N_1  (Ctj*  -  1)  +  1  (85) 

u 

s 
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where  C  is  a  coefficient  to  be  determined.  It  is  clear 
that  Eq.  (85)  is  a  polynomial  of  degree  N  and  that  the 
outer  edge  condition  u/us(t)*  =  l)  =  1  is  automatically 
satisfied.  The  non  slip  condition  at  r*  =  0  is  also 
automatically  satisfied.  It  may  be  shown  by  differentiation 
that  the  N  -  3  outer  edge  conditions  given  by  Eq.  (80)  are 
also  automatically  satisfied.  The  remaining  boundary 
condition  is  the  second  condition  of  Eq.  (79)  which  can 
be  used  to  determine  C  with  the  result 

c  _  x*  -  (H  -  1)(M  -  2)  (86) 

2(H  -  1> 


Substituting  Eq.  (86)  into  (85)  gives 


— —  =  (1  “ 
u 


m 


X*  -  (N-l ) (N-2 ) 
2 (N  -  1) 


-  1 


+  1 


(87) 


as  an  equivalent  form  of  the  Nth  degree  polynomial  profile 
which  satisfies  all  the  required  boundary  conditions.  Ex¬ 
panding  the  term  (1  -  q*)N-l  by  the  binomial  theorem  then 
gives  a  polynomial  of  degree  N.  Equating  the  coefficients 
of  successive  powers  of  r|*  from  Eqs.  (78)  and  (87)  gives 
exact  recurrence  relations  for  the  an  with  the  result 


an  -  <-!> 


n+1 


N  X*  ~  (N-l)  (N-2)  _  fN-1| 
n_1/  2  (N  -  1)  l  N  > 


for  n  =  1,  2,  3  ...N-l  and 


aN  = 


(-1) 


N+1  f  X*  -  (N-l)  (N-2 
l  2 (N  -  1) 


)' 


where 


N-l 

n-l 


(N-l)!  n 
nj  (N-n); 


N-ll  =  (N-l);  (N-n) 
n  nl  (N  -  n) ; 


(88) 


(89) 
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The  wall  shear  stress,  appearing  in  Eq.  (73)  can 

now  be  explicitly  related  to  a*  and  N  with  the  result 


T  =  |X  — — 

W  dz 


z=0 


(- 


2 (N-l) 


+  S) 


(90) 


Physically  meaningful  boundary  layer  results  can  be 
obtained  for  A*  values  in  a  well  defined  range.  For 
sufficiently  small  values  of  A*,  separation  occurs  where¬ 
as  for  large  values  of  A*,  velocity  overshoot  in  the 
boundary  layer  occurs.  The  condition  for  boundary  layer 
separation  at  the  wall  is 


=  0 

z=0 


which  gives 


=  -  N(N  -  1) 


(91) 


Velocity  overshoot  occurs  in  the  boundary  layer  if  —  =  0 
at  some  intermediate  value  of  z  between  0  and  A.  dz 
Applying  this  condition  gives 


Xmax  *  N<N  "  D 


(92) 


Thus  the  physically  allowable,  range  for  A*  becomes 


-  N(N  -  1)  <  A*  <  N(N  -  1) 


(93) 


where  the  negative  lower  limit  corresponds  to  separation 
and  the  positive  upper  limit  corresponds  to  velocity  over¬ 
shoot  . 


With  the  expressions  for  an  now  known  in  explicit  form, 
the  displacement  and  momentum  thicknesses  may  be  evaluated 
from  Eqs.  (74)  and  (75)  with  the  result 


6*=  3N(N  -  1)  -  A* 
A  2N(N2  -  1) 


28 


AEDC-TR- 67-79 


S*=  (W-1) _ (N-2)C _ c2 

A  N(2N-1)  N(N+1) (2N-1)  N(2N+1 ) (2N-1 ) 


(95) 


where 

c  =  X*  -  (N-1HH-2) 
2 ( N- 1 ) 


Combining  Eqs.  (73),  (77),  (94)  and  (95)  leads  to  a 
single  differential  equation  involving  X*  and  0*  which 
can  be  written  in  the  form 


dX* 

d@* 


f(x* ,  e*) 


(96) 


The  derivation  of  Eq.  (96)  is  quite  lengthy  and  will  not 
be  presented  here,  Eq.  (96)  is  non  linear  and  must  be 
integrated  numerically.  Once  the  numerical  relationship 
between  X*  and  9 *  is  established,  one  may  then  obtain 
numerically  the  variations  of  &*,  X*,  5*,  tw,  us>  and  P 
with  x.  Of  particular  interest  here  is  the  variation 
of  us,  the  free  stream  velocity,  which  gives  a  measure  of 
the  entrance  length.  The  influence  of  the  degree  of  the 
polynomial,  N,  on  the  entrance  length  is  shown  in  Table  2. 
The  entrance  length  is  defined  as  the  distance  at  which 
the  centerline  velocity  becomes  95%  of  the  fully  developed 
centerline  velocity.  From  the  table  it  is  clear  that 
increasing  N  from  4  to  60  does  not  have  a  strong  effect 
on  the  calculated  entrance  length.  The  influence  of  N  is 
more  significant  for  smaller  values  of  Ha,  and  for  Ha  =  1, 
the  difference  in  the  calculated  entrance  lengths  for  N  =  4 
and  N  =  60  is  8%.  The  entrance  length  decreases  slightly 
as  N  increases. 

3.11  THE  INFLUENCE  OF  VARIABLE  CONDUCTIVITY  ON  MHD  CHANNEL  FLOW 


The  analysis  presented  thus  far  has  been  based  on  the 
assumption  of  constant  fluid  properties,  including  the 
electrical  conductivity.  The  electrical  conductivity  is  a 
strong  function  of  temperature,  however,  and  any  flow  with 
large  temperature  variations  will  have  a  spatially  non 
uniform  conductivity. 

The  direct  solution  of  variable  conductivity  MHD  flows 
requires  the  simultaneous  solution  of  the  energy  equation 
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and  the  equation  of  motion.  An  iterative  procedure  is 
normally  used  in  which  a  tentative  solution  to  the  equation 
of  motion  is  obtained  using  an  assumed  spatial  variation  of 
the  conductivity.  Using  this  solution  in  the  energy  equation 
then  gives  a  first  order  solution  for  the  temperature  pro¬ 
file.  From  the  temperature  profile  and  a  prescribed  temper¬ 
ature-conductivity  variation,  a  corrected  conductivity 
profile  can  be  obtained.  Using  the  corrected  conductivity 
profile  in  the  equation  of  motion  allows  the  procedure  to 
be  repeated  as  many  times  as  desired. 

The  effectiveness  of  the  iteration  procedure  depends 
on  the  difficulty  involved  in  solving  the  equation  of 
motion  with  a  spatial  variation  in  conductivity.  Exact 
solutions  can  be  obtained  for  only  a  limited  number  of 
conductivity  variations  for  fully  developed,  incompressible 
flow.  Kieffer  (Ref.  20)  obtained  exact  solutions  for  three 
cases  of  conductivity  variation,  namely  a  linear  variation, 
a  parabolic  variation,  and  a  step  function  variation. 

These  solutions  appear  to  be  the  only  closed  form  solutions 
in  the  literature  for  variable  conductivity  channel  flow. 

In  the  present  investigation,  two  techniques  of  analysis 
of  variable  conductivity  channel  flows  are  presented.  Both 
of  these  methods  are  capable  of  handling  an  arbitrary  spatial 
variation  of  conductivity  in  the  z-direction.  The  two 
methods  employed  are  the  use  of  Fredholm  integral  equation 
theory  and  the  B.  G.  Galerkin  method.  The  analysis  is 
restricted  to  the  fully  developed  flow  region. 

The  Application  of  Fredholm  Integral  Equation  Theory  to 
Variable  Conductivity  MHD  Channel  Flows 

The  equation  of  motion  for  fully  developed  flow  in 
the  absence  of  Hall  current  may  be  written  from  Eq.  (5)  as 


0  = 


- 

dx 


+  “  X?  + 


UBZ> 


(97) 


The  origin  z  =  0  as  shown  in  Fig.  la  is  located  at  the 
centerline  of  the  channel.  For  the  present  analysis  it 
is  more  convenient  to  locate  the  origin  at  the  lower  sur¬ 
face.  Also,  non  dimensionalizing  the  velocity  in  terms 
of  the  pressure  gradient  is  convenient.  Thus  the  following 
dimensionless  quantities  are  defined. 


(2h)2  d£ 
ji  dx 
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r 


a 


Ha 


(2h)2o  B  ‘ 

*2  _  _ O  Z 


c  B  E 

«  -  °.  2-2.  n*  =  — 

y  d£  2h 

dx 


08) 


where  aQ  is  a  constant  reference  value  of  conductivity. 

The  notation  ap  indicates  the  velocity  profile  obtained  by 
using  the  Fredholm  integral  equation  theory.  The  dimen¬ 
sionless  form  of  Eq.  (97)  becomes 


d2a„ 

- *  -  Ha*2r(r1*)aT?  +  4>*r(q*)  +1  =  0  (99) 

dq*2  F  y 


Eq.  (99)  can  be  converted  to  an  integral  equation  by  in¬ 
tegrating  the  equation  twice,  subject  to  the  boundary 
conditions  aF(q*  =  0)  =  0  and  aj^q*  =  .l)  =  0,  The  resulting 
integral  equation  may  be  written 


1 

ap(q*)  =  f  (q* )  +  Ha*2  /  K(q* , t )a* (t)dt 

0 

where 

f(q*)  =  /  ( t  -  r)*)[<fc*r(t)  +  l]dt 

0 

1 

+  q*  /  (1  -  t)[<D*r(t)  +  l]dt 
0 


K(q*,t)  =  t(q*  -  l)T(t)  0  <  t  £  q* 


K(q*,t)  =  q*(t  -  l)r(t)  q*  <  t  <  1 


(100) 


(101) 


•  (102) 
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Eq.  (100)  is  a  Fredholm  integral  equation  with  kernel 
K(n*,t)  defined  by  Eq.  (102).  To  arrive  at  the  standard 
form  of  the  Fredholm  integral  equation,  it  is  necessary 
to  define  the  kernel  differently  in  different  regions  of 
the  interval  as  shown  by  Eq.  (102).  The  kernel  is  con¬ 
tinuous,  however,  over  the  entire  interval. 


The  fact  that  the  kernel  vanishes  at  the  end  points 
of  the  interval  t  =  0  and  t  =  1  suggests  a  sine  series 
expansion  of  the  form 


K(n*,t) 


oo 

^  an(r]*)sin(nirt)dt 
n=l 


(103) 


The  functions  an(r|*)  may  be  evaluated  using  the  ortho¬ 
gonality  of  the  sine  function  with  the  result 

1 

an(T|*)  =  2  /  K(t;*  ,  t )sin (niTt)dt  (104) 

0 

i 

* 

Substituting  the  definition  of  the  kernel  from  Eq.  (102). 
into  Eq.  (104)  gives 


an(rj*)  =  2 [  /  t(q*  -  l)r(t)sin(nirt)dt 
0 

1 

+  /  q*(t  -  l)r(t)sin(mrt)dt]  (105) 

n* 


For  a  degenerate  kernel  consisting  of  a  finite  number 
of  terms,  an  exact  solution  to  the  Fredholm  integral 
equations  may  be  obtained.  A  degenerate  kernel  is  defined 
as  a  kernel  in  which  each  term  consists  of  a  product  of 
two  functions,  each  function  depending  on  only  one  of  the 
variables  in  the  kernel.  K(rj*,t)  as  given  by  Eq.  (103) 
is  seen  to  be  degenerate,  and  if  the  Fourier  series  is 
truncated  after  a  finite  number  of  terms,  an  exact  solution 
to1 the  integral  equation  can  be  found.  In  the  present 
analysis  the  series  given  by  Eq,  (103)  is  truncated  after 
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M  terms.  This  results  in  an  integral  equation  which  is 
an  approximation  to  the  original  integral  equation.  How¬ 
ever,  an  exact  solution  to  the  approximate  integral  equation 
can  be  obtained. 


The  solution  details  for  the  Fredholm  integral  equation 
are  described  in  Lovitt  (Ref.  21)  and  only  the  results 
will  be  presented  here.  With  the  kernel  of  Eq.  (103)  trun¬ 
cated  after  M  terms,  the  solution  to  Eq.  (100)  is 


ctj^r)*) 


M 


f (r* )  +  Ha 


*2 


I 

n=l 


C*a 
n  n 


(h*) 


(106) 


The  C*  coefficients  are  the  solutions  to  the  set  of  linear 
algebraic  equations 


M 


Ha 


+  2 


E 


i=l 


B.  C. 
in  l 


In  n  -  1 , 2 , . . . M 


(107) 


I 


n 


1  +  (-D 

n37r3 


n+1 


2  2 
n  V 


1 

/  r(t)sin(nirt)dt 
0 


(108) 


1 

B.  =  -  — — —  /  T  (t )sin (irt) sin (nirt  )dt  (109) 

in  22 


Prescribing  a  conductivity  profile  T(t)  permits 
evaluation  of  the  terms  In  and  B-^n.  It  should  be  noted 
that  T(t)  appears  only  within  integrals  in  the  above  ex¬ 
pressions  and  thus  a  conductivity  profile  in  tabular  form 
could  also  be  handled. 

Numerical  results  were  obtained  for  three  cases  of 
conductivity  profile,  uniform,  linear,  and  parabolic.  The 
assumed  conductivity  variations  for  the  three  cases  are 

Case  I  r(r)*)  =  1 


33 


AEOC-T  R -67-79 


Case  II 


Case  III 


r(Ti*)  =  2r|*  0  <  q*  £  0.5 

r  (tj*)  -  2(1  -  if)  0.5  <  q*  <  1.0 


rCn*)  -  4<-n*  -  r,*2) 


In  cases  II  and  III,  the  reference  conductivity  oQ  has 
been  taken  as  the  centerline  value.  These  profiles  were 
chosen  because  exact  solutions  to  the  differential  equation 
of  motion  can  be  obtained  with  which  to  compare  the  ap¬ 
proximate  results  of  the  present  analysis. 

Velocity  profiles  for  the  three  cases  of  conductivity 
variation  are  shown  in  Figs.  22  through  29.  Each  of  these 
curves  is  based  on  a  truncation  of  the  kernel  expansion 
given  by  Eq.  (103)  after  three  terms,  i.e.  M  =  3.  Shown 
also  are  the  results  from  an  exact  solution  of  the  dif¬ 
ferential  equation.  The  agreement  between  the  exact  solution 
to  the  differential  equation  and  the  solution  to  the  integral 
equation  using  a  kernel  series  truncated  after  three  terms 
is  excellent. 

To  investigate  the  effect  of  the  number  of  terms  re¬ 
tained  in  the  kernel  expansion,  M  was  varied  from  3  to  20 
with  several  intermediate  values.  In  all  cases,  the  dis¬ 
crepancy  between  the  3  and  20  term  solutions  never  ex¬ 
ceeded  3%  and  in  most  cases  was  less  than  1%. 

The  application  of  the  Fredholm  integral  equation 
theory  to  variable  conductivity  MHD  channel  flows  appears 
to  be  a  powerful  technique  for  obtaining  an  accurate 
approximation  to  the  fully  developed  velocity  profile  for 
an  arbitrary  spatial  variation  of  conductivity.  Since  the 
conductivity  profile  T(q*)  always  appears  in  an  integral, 
tabulated  conductivity  data  could  also  be  used,  requiring 
only  that  the  integrals  be  evaluated  numerically. 

The  Application  of  the  B.  G.  Galerkin  Method  to  Variable 
Conductivity  MHD  Channel  Flows 

Another  method  suitable  for  analyzing  variable  con¬ 
ductivity  channel  flows  is  the  B.  G.  Galerkin  method.  The 
B.  G.  Galerkin  method  is  discussed  extensively  in  (Ref.  22) 
and  only  the  essentials  of  the  method  as  it  applies  to  the 
present  problem  will  be  discussed  here.  It  is  convenient 
to  non  dimensionalize  Eq.  (97)  in  a  manner  slightly  different 
from  that  used  in  the  application  of  the  Fredholm  integral 
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equation  theory.  Hie  following  dimensionless  quantities 
are  defined. 


|x  dx 


h2a  B  2 

2  O  Z 


Ha  = 


a  BE 
o  z  y 

dP 

dx 


z 

T)  =  “ 
h 


(111) 


The  origin  z  =  0  will  once  again  be  taken  at  the  center- 
line  of  the  channel.  Since  the  B.  G,  Galerkin  method  is 
most  amenable  to  two  point  boundary  conditions,  it  is 
convenient  to  obtain  the  solution  over  the  entire  channel 
instead  of  the  half  channel.  Thus  the  range  of  q  will  be 
-  1  ^  T)  <  1 .  The  dimensionless  form  of  Eq.  (97)  becomes 


d2ar 

- ±  -  Ha2T (q)ar  =  -  [<&*r(r})  +  1]  (112) 

dr)2  G  y 


with  boundary  conditions  0^(1)  -  -1)  =  0  and  aG(q  =  1)  =  0. 

The  essential  idea  of  the  B,  G.  Galerkin  method  is  to 
construct  an  approximate  solution  to  Eq.  (112)  as  a  sum  of 
functions  each  of  which  satisfies  the  boundary  conditions. 
Let  Gk(q)  be  a  set  of  continuous  functions  which  satisfy 
the  boundary  conditions 


Gk(T)  =  -  1)  =  Gk(q  =  1)  =  0  k  =  1,2,... S  (113) 


An  approximate  solution  to  Eq.  (112)  can  be  constructed 
from  a  linear  combination  of  the  Gk(~})  functions  as 
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s 

°G  =  ]T  bnGn(T|)  (114) 

n=0 


From  Eq.  (113)  it  is  clear  that  the  boundary  conditions 
of  Eq.  (112)  are  satisfied  by  the  solution  given  by  Eq. 
(114).  It  remains  to  determine  the  coefficients  b^  such 
that  the  differential  equation  is  satisfied.  Substituting 
Eq.  (114)  into  Eq.  (112)  gives 

S 

Z 

n=0 


Multiplying  Eq.  (115)  by  each  function  of  the  set  and 
integrating  over  the  range  of  q  gives 

1 

bn  f  rGkGndri 
-1 


k=0,l,2,...S  (116) 

1 


=  -<!>*/  r-Gkdn  -  /  Gkdn 


b  G 
n  n 


-  Ha2r 


z 

n=0 


b  G 
n  n 


-  cv  + 


i] 


(115) 


Once  a  set  of  approximating  functions  Gd(t])  is  chosen, 
the  definite  integrals  in  Eq.  (116)  can  be  evaluated.  Eq. 

(116)  represents  a  set  of  S  4-  1  linear  equations  for  the 
S  +  1  unknown  coefficients  bn.  With  the  bn  determined,  the 
approximate  solution  to  the  differential  equation  is  completed. 

For  the  case  of  fully  developed  MHD  channel  flow,  a  set 
of  approximating  functions  which  satisfy  the  boundary  con¬ 
ditions  may  be  taken  as 


Gn  =  (1  -  Ti2)T}n  n=0, 1, 2. .  . S  (117) 

The  approximate  velocity  profile  for  this  set  of  functions 
becomes 


36 


AEDC-T  R -67*79 


s 


n=0 


(118) 


Substituting  Eq.  (117)  into  Eq „  (116)  and  performing 

the  integrations  gives  the  following  system  of  linear 
equations  for  the  bn’s- 

b  [ [1  -  (-l)k+1l  ,  Ha2 
°  \  (k+1 ) (k+3 )  4 


1(0, k)  + 


b 

i 


[1  +  (-l)k+1] 

(k+2) (k+4) 


+  Sal  i(i,k)| 


+ 


E" 


[i  +  (-nn+kj 


n(k+l)  +  k(n+l)  -  1 
(n+k-1) (n+k+1) (n+k+3) 


Ha2  \\ 

—  I  (n,k) 

4 

// 


1  fl  -  (-1)K+J-1  'V 

ill - - L - I(k)  k=0, 1,2.  .  .  S  (119) 

2  (k+1) (k+3)  4 


where  the  integrals  I(n,k)  and  I(k)  are  given  by 

1 

I  (n,k)  =  /  r (q)  (r,2  -  l)2qn+kdr, 

-1 

1 

I(k)  =  /  F(t))(t]2  -  l)r^kdr| 

-1 


(120) 


For  any  given  F  (rj)  conductivity  profile  the  integrals 

of  Eqs.  (120)  can  be  evaluated  either  analytically  or 

numerically  and  the  linear  set  of  equations  given  by  Eq. 

(119)  solved  for  the  unknown  constants  b  . 

n 
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Velocity  profiles  for  Cases  I  and  III  were  also  ob- 
tained  by  the  B.  G.  Galerkin  method  using  a  set  of  20 
approximating  functions  as  defined  by  Eq.  (117)  with  S  =  19. 
The  agreement  between  this  approximate  solution  and  the 
exact  solution  to  the  differential  equation  was  excellent. 
For  Case  I,  the  maximum  difference  between  the  exact  and 
approximate  solutions  was  0.003%  whereas  for  Case  II,  the 
maximum  difference  was  less  than  1%. 

In  addition  another  non  uniform  conductivity  profile 
was  investigated  by  the  B.  G.  Galerkin  method.  Because  of 
the  strong  dependence  of  conductivity  on  temperature,  in 
the  vicinity  of  a  highly  cooled  wall  the  conductivity  may 
drop  to  essentially  zero.  A  profile  which  exhibits  a  zero 
conductivity  layer  on  each  wall  is  given  by  the  expression 


Case  IV 

r(q)  =  0 

V| 

rH 

1 

-ZJ 

'w' 

II 

“  i  < 

r(r,)  =  b 

i  < 

n  <  i 
<  I 

T)  ^  1  (121) 


The  profile  given  by  Eq.  (121)  is  discontinuous  with  a 
central  portion  of  uniform  conductivity  of  thickness  2£ 
and  zero  conductivity  layers  of  thickness  £  on  each  B  wall. 

Numerical  results  were  obtained  for  Case  IV  for  £  - 
0.5  which  corresponds  to  a  zero  conductivity  layer  of 
thickness  one  quarter  of  the  channel  height.  The  resulting 
velocity  profiles  are  shown  in  Figs.  30  and  31  for  short 
circuit  and  open  circuit  conditions  and  for  Hartmann  numbers 
Ha  =  6  and  10.  Shown  also  for  comparison  are  the  profiles 
corresponding  to  uniform  conductivity. 

The  profiles  shown  in  Fig.  30  are  for  the  short  cir¬ 
cuit  case,  i.e.  Ey  =  0  =  $*.  A  large  velocity  overshoot 
occurs  in  which  the  maximum  velocity  is  not  on  the  center- 
line  of  the  channel  but  occurs  near  the  wall.  The  profiles 
shown  in  Fig.  31  are  for  the  open  circuit  condition.  The 
velocity  overshoot  does  not  occur  for  this  loading  con¬ 
dition,  but  the  profiles  are  qualitatively  quite  different 
from  the  uniform  conductivity  profiles.  In  particular  the 
shear  stress  at  the  wall  is  smaller  for  the  non  uniform 
conductivity  case.  Also  the  width  of  the  central  core  of 
essentially  uniform  velocity  is  decreased  due  to  the  non 
uniform  conductivity. 
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It  should  be  noted  that  the  velocities  in  Figs.  30  and 
31  are  normalized  in  terms  of  the  pressure  gradient  instead 
of  the  mean  velocity.  This  means  that  for  a  given  value 
of  Ha,  the  curves  for  uniform  and  non  uniform  conductivity 
correspond  to  different  values  of  the  mean  velocity  but  the 
same  value  of  the  pressure  gradient.  The  areas  under  the 
curves,  which  would  give  the  mean  velocity,  are  therefore 
different  at  the  same  value  of  Ha.  A  relation  between  the 
mean  velocity  and  the  pressure  gradient  is  derived  in  the 
next  section. 

Operational  Characteristics  of  the  Channel 

In  analyzing  variable  conductivity  channel  flows,  it 
is  convenient  from  a  mathematical  viewpoint  to  non 
dimensionalize  the  velocity  in  terms  of  the  pressure  gra¬ 
dient  instead  'of  the  mean  velocity.  This  results  in  a 
loading  parameter  $y  defined  by  Eqs.  (98)  and  (111)  in  terms 
of  the  pressure  gradient.  For  uniform  conductivity,  the 
loading  parameter  is_normally  defined  in  terms  of  the  mean 
velocity  as  =  Ey/UBZ.  It  is  desirable  to  obtain  a  re¬ 
lation  between  <£*  and  $y  , 

The  mean  velocity  is  defined  as 

h 

U  =  —  /  udz  (122) 

2h-h 

which  may  be  rewritten  in  terms  of  dimensionless  quantities 
as 


U  = 


1 


h2  dP  f  , 
- J  ad-] 

2\i  dx  n 


(123) 


Eq.  (123)  relates  the  mean  velocity  to  the  pressure  gradient 
in  terms  of  the  integrated  velocity  profile,  the  velocity 
being  non  dimensionalized  in  terms  of  the  pressure  gradient. 
For  specified  values  of  Ey  and  Bz,  the  relation  between 

and  <t*  becomes  y 

y 


y  2  y  , 


(124) 


where  Ha2  is  the  Hartmann  number  defined  by  Eqs,  (111). 
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The  other  physical  quantities  of  interest  are  the  open 
circuit  values  of  <I>  and  <t>y .  For  uniform  conductivity , 
the  open  circuit  value  of  is  1  but  will  be  different 
from  1  for  non  uniform  conductivity.  The  open  circuit 
condition  can  be  established  by  considering  the  total  cur¬ 
rent  flowing  in  the  y  direction.  The  total  y  component  of 
current  per  unit  length  of  the  channel  in  the  x  direction  is 

h  1 

Iy  -  J  Jydz  =  h  /  a (Ey  -  uBz)dq  {125) 

-h  -1 

Writing  Eq.  (125)  in  terms  of  dimensionless  quantities 
gives 

1 

I  =  -  "  —  /  r(3>*  -  Ha2a)dTj  (126) 

y  B  dx  _  y 


The  open  circuit  condition  corresponds  to  Iy  =  0. 
for  with  I  =0  gives 

y  y  1 

Ha2  /  Tadtj 


/  ran 

-1 


Solving 


(127) 


A  value  of  <£*  must  be  assumed  to  obtain  a  solution 
for  a  by  either  the  Fredholm  integral  equation  method  or 
the  B.  G.  Galerkin  method.  Since  the  right  side  of  Eq. 
(127)  depends  on  the  value  of  <f*  assumed,  the  equation 
must  be  solved  by  trial  and  error  for  $y)QC>  That  is, 
for  a  fixed  value  of  Ha,  the  right  side  or  Eq.  (127)  is 
evaluated  for  different  assumed  values  of  4>y.  The 
particular  value  of  $y  for  which  the  right  side  of  Eq. 
(127)  equals  <3>y  is  the  open  circuit  value  |ty)oc*  Com¬ 
bining  Eqs.  (126)  and  (127)  gives^ 

/  rad-rj 


oc 


-1 

2  - 

1  1 

/  adq  /  Tdr| 

-1  -1 


(128) 
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It  is  clear  from  Eq.  (128)  that  <±>  )oc  =  for  r  =  1 , 
corresponding  to  uniform  conductivity. 

Values  of  <£y)oc  calculated  from  Eq.  (128)  are  tab¬ 
ulated  in  Table  3  for  several  values  of  Ha  for  Case  IV. 

The  value  of  $  )  is  seen  to  decrease  slightly  with  in- 
TT  y  oc 
creasing  Ha.  J 

The  relations  contained  in  Eqs.  (122)  through  (128) 
are  valid  for  either  the  Fredholm  integral  equation  theory 
or  the  B.  G.  Galerkin  method  of  analysing  variable  con¬ 
ductivity,  fully  developed,  MHD  channel  flow. 
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IV 

EXPERIMENTAL  INVESTIGATION 


4.1  EXPERIMENTAL  FACILITY 

There  is  very  little  experimental  data  available  in 
the  open  literature  for  the  entrance  region  of  channels, 
even  for  non  MHD  flows.  No  experimental  velocity  profile 
data  for  MHD  entrance  region  flows  appear  to  be  available. 

In  order  to  obtain  such  data  for  determining  the  validity 
of  the  analytical  results  presented  in  previous  sections, 
a  new  experimental  MHD  facility  was  designed  and  constructed 
as  part  of  the  present  investigation.  The  design  and  fab¬ 
rication  of  the  facility  was  done  by  the  Engineering  Support 
Facility  of  ARO,  Inc.  The  experiments  were  conducted  by  the 
Experimental  Research  Group,  Technical  Staff,  OMD  of  ARO, 

Inc  . 


A  schematic  drawing  of  the  flow  loop  is  shown  in  Fig. 

32.  The  loop  is  closed  with  a  pump  used  for  circulating 
the  working  fluid.  The  major  components  of  the  loop  as 
shown  in  Fig,  32  consist  of  a  test  section,  a  downstream 
section  joining  the  test  section  to  the  piping,  flowmeters, 
a  circulating  pump,  a  pump  bypass  loop  for  flow  control,  a 
transition  section  from  the  piping  to  the  plenum  chamber, 
a  plenum  chamber,  and  appropriate  valving  for  controlling 
the  flow  rate.  The  maximum  capacity  of  the  circulating 
pump  is  approximately  40  gpm  which  gives  a  maximum  test 
section  Reynolds  number  of  11,600.  The  channel  and  test 
section  are  designed  for  use  with  an  electrolytic  solution 
of  potassium  chloride  (KC1)  as  the  working  fluid.  The  test 
section  is  made  of  plexiglass  and  the  circulating  piping  is 
stainless  steel.  The  system  has  venting  ports  and  valves 
at  several  locations  on  the  top  side  of  the  loop,  particularly 
in  the  vicinity  of  abrupt  cross  sectional  area  changes,  for 
removing  entrapped  gas  from  the  system. 

A  cross  sectional  view  of  the  test  section  in  the  x-z 
plane  is  shown  in  Fig.  33.  The  test  section  is  rectangular 
and  is  constructed  as  a  channel  within  a  channel.  The  inner 
channel  walls  in  the  x-y  plane,  the  so  called  B-walls  since 
they  arc  perpendicular  to  the  magnetic  field,  are  constructed 
with  a  sharp  leading  edge  to  provide  a  well  defined  inlet 
section  with  a  reasonably  flat  initial  velocity  profile. 

The  inner  channel  is  provided  with  a  total  of  22  static 
pressure  measuring  stations  one  inch  apart,  the  first  static 
pressure  station  being  located  2  inches  from  the  inlet.  At 
each  pressure  measuring  station  in  the  channel,  the  pressure 
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taps  on  opposite  sides  of  the  channel  are  interconnected 
with  a  common  pressure  line  leading  to  the  manometer. 

The  dimensions  of  the  inside  channel  are  3/4  x  5  1/2 
inches  giving  an  aspect  ratio  of  7  1/3.  This  value  of  aspect 
ratio  should  be  large  enough  to  closely  approximate  a  parallel 
plate  geometry  in  the  center  portion  of  the  channel.  The 
two  surfaces  of  the  inner  channel  in  the  x-z  plane  are 
electrode  surfaces.  The  electrodes  are  continuous  platinum 
sheets  0.001  inch  thick  which  are  cemented  to  the  channel 
walls.  Electrical  leads  are  attached  to  both  electrode 
surfaces . 

The  pressure  lines  from  the  22  static  pressure  measuring 
stations  are . individually  separated  from  a  common  manifold 
by  a  needle  valve.  The  pressure  from  any  station  can  be 
communicated  to  the  manifold  by  opening  the  appropriate 
valve,  the  remaining  valves  being  closed.  A  pressure  line 
from  the  manifold  is  communicated  to  one  leg  of  a  micro¬ 
manometer  . 

A  pitot  probe  with  traversing  mechanism  is  installed 
for  measuring  velocity  profiles.  The  tip  of  the  probe  is 
located  at  the  last  static  pressure  measuring  station  which 
is  23  inches  from  the  entrance  to  the  inner  channel.  The 
probe  can  be  traversed  across  the  channel  in  the  z  direction 
and  can  be  positioned  to  an  accuracy  of  0.001  inch  by  means 
of  a  micrometer  screw.  The  probe  cannot  be  moved  in  the  x- 
direction.  By  adjusting  the  flow  rate  and  thereby  the 
channel  Reynolds  number,  the  length  of  the  entrance  region 
can  be  varied.  In  this  way  the  pitot  tube  can  sense  different 
portions  of  the  entrance  region  flow  by  remaining  at  the 
same  x-position. 

The  magnetic  field  is  provided  by  an  electromagnet 
with  rectangular  pole  pieces.  The  pole  pieces  are  6  x  30 
inches  with  several  sets  of  different  thickness  available 
to  vary  the  gap  width.  The  maximum  gap  width  obtainable  is 
2  inches.  The  magnet  dimensions  are  such  that  the  entire 
test  section  can  be  inserted  in  the  air  gap.  The  magnet  is 
powered  by  D.C,  electric  welders.  The  maximum  power  rating 
of  the  magnet  is  600  K.W.  At  maximum  power,  a  field  strength 
of  22,000  gauss  can  be  produced  in  a  2  1/2  inch  air  gap. 
Magnetic  field  uniformity  has  been  measured  to  indicate  less 
than  5%  non  uniformity  along  the  centerline  of  the  magnet 
over  the  central  24  inches  of  pole  piece  length.  The  test 
section  inlet  is  located  4  inches  from  the  magnet  so  that 
the  entire  entrance  region  of  the  channel  is  in  a  magnetic 
field  with  less  than  5%  variation  in  the  field  strength 
uniformity . 
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A  photograph  of  the  test  section,  flow  loop,  and 
micromanometer  is  shown  in  Fig.  34.  In  the  interest  of 
clarity,  the  electromagnet  is  not  shown. 

4.2  INSTRUMENTATION 

The  principal  measurements  involved  in  the  experiment 
are  the  fluid  flow  rate,  the  static  pressure  distribution  in 
the  channel,  and  the  velocity  profile  across  the  3/4  inch 
dimension  of  the  inner  channel  in  the  z-direction,  i.e.  the 
velocity  profile  in  a  plane  normal  to  the  B-wall.  The  flow 
rate  is  measured  by  variable  area  flow  meters,  calibrated 
by  the  Instrumentation  Branch  of  Engineering  Support  Facility 
of  ARO,  Inc.  Two  flow  meters  are  used,  one  for  low  flow 
rates  and  the  other  for  high  flow  rates.  The  pitot  probe, 
located  at  the  last  static  pressure  station  23  inches  from 
the  channel  inlet,  is  used  to  measure  velocity  profiles  by 
communicating  the  probe  pressure  to  one  leg  of  a  micro¬ 
manometer  and  the  channel  wall  static  pressure  to  the  other 
leg.  In  the  vicinity  of  the  channel  wall,  the  velocity 
head  to  be  measured  is  extremely  small  as  will  be  demon¬ 
strated  subsequently. 

To  maintain  laminar  flow  in  the  channel,  the  Reynolds 
number  must  be  kept  sufficiently  low.  Using  a  liquid  with 
the  density  of  water  as  the  working  fluid,  the  relationship 
between  the  velocity  head  and  the  Reynolds  number  may  be 
written  as 


H  = 


PUg 

2 


1.92  x  10  e(Re)g  (inches  of  water) 


(129) 


The  characteristic  dimension  used  in  the  definition  of  the 
Reynolds  number  is  the  channel  half  height  in  the  z-direction, 
h  =  0.375  inches.  Eq,  (129)  can  be  used  to  calculate  the 
velocity  head,  i.e.  difference  between  total  and  static 
pressure,  to  be  measured.  For  example  a  value  of  Re  -  5,000 
gives  H  =  0.48  inch  of  water  corresponding  to  the  mean  flow 
velocity  in  the  channel.  As  the  pitot  probe  approaches  the 
channel  wall,  the  velocity  head  decreases  to  zero  at  the 
wall.  The  minimum  velocity  measurable  would  be  with  the 
pitot  probe  touching  the  channel  wall.  The  outside  diameter 
of  the  pitot  probe  is  0.063  inch  which  means  the  closest 
position  at  which  a  velocity  measurement  can  be  made  is 
0.031  inch  from  the  wall.  For  fully  developed  non  MHD  flow, 
the  velocity  head  at  this  location  in  terms  of  the  Reynolds 
number  is 

Hnin  -  1.11  x  10-9(Re)2  (inches  of  water)  (130) 
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For  Re  =  5,000,  Hmin  =  0.028  inch  of  water.  Thus  in  the 
vicinity  of  the  wall,  extremely  small  pressure  measurements 
are  involved  in  determining  the  velocity  head. 

Initial  efforts  were  made  to  measure  the  velocity  head 
using  a  commercially  available  Hass  micromanometer.  This 
instrument  is  capable  of  measuring  a  pressure  difference  as 
small  as  0.0001  inch  of  mercury.  The  Hass  micromanometer 
is  in  principle  a  very  sensitive  U-tube  manometer.  Two 
glass  cylinders  containing  the  manometer  fluid  are  inter¬ 
connected  through  a  tube.  Applying  a  pressure  differential 
to  the  free  surfaces  of  the  manometer  fluid  in  each  cylinder 
results  in  a  deflection.  The  fluid  level  in  each  cylinder 
is  determined  by  rods,  sharpened  to  a  needle  point,  actuated 
in  the  vertical  direction  by  micrometer  screws.  When  the 
needle  points  touch  the  liquid  surface,  a  small  dimple  is 
formed  on  the  surface  which  can  be  observed  through  a  mag¬ 
nifying  glass.  With  no  pressure  difference  applied  across 
the  manometer,  the  micrometers  are  read  for  the  rod  positions 
at  which  the  needle  points  touch  the  surfaces.  When  a  pressure 
differential  is  applied,  the  liquid  surfaces  deflect  and 
the  deflections  of  the  two  surfaces  are  determined  by  the 
rod  travel  required  for  the  needle  points  to  touch  the 
surfaces  again. 

Considerable  difficulty  was  experienced  in  obtaining 
reproducible  pressure  measurements  with  the  micromanometer. 

The  difficulty  was  attributable  to  leaks  in  the  manometer 
and  connecting  lines  and  to  the  extremely  slow  response  time 
of  the  manometer.  Generally  four  to  five  hours  were  required 
to  obtain  one  data  point  due  to  the  slow  response  of  the 
manometer. 

Several  weeks  time  were  devoted  to  attempting  to 
solve  the  leak  problem.  After  ascertaining  that  all  valves 
and  connections  on  the  pressure  lines  were  leak  free,  it  was 
finally  discovered  that  the  tygon  tubing  was  defective,  being 
sufficiently  porous  to  allow  leaks  to  occur  while  waiting 
for  the  manometer  to  stabilize.  Consequently,  all  tygon 
tubing  was  replaced  by  metal  tubing. 

Leaks  were  also  detected  at  the  "O"  ring  seals  on  the 
movable  rods  in  the  micromanometer.  New  ,T0'T  rings  of  smaller 
diameter  were  installed  and  the  packing  nuts  were  tightened 
as  tight  as  possible  and  still  permit  movement  of  the  rods. 

By  retightening  the  packing  nuts  before  each  run,  it  was 
possible  to  obtain  a  limited  number  of  pressure  measurements. 

Because  of  the  severe  problems  experienced  in  using 
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the  micromanometer  system,  an  alternate  method  of  pressure 
measurement  was  explored.  A  considerable  effort  was  ex¬ 
pended  in  attempting  to  use  an  electrokinetic  pressure 
transducer.  The  theory  of  this  device  is  discussed  in 
(Ref.  23).  The  electrokinetic  transducer  has  the  advantage 
of  an  extremely  fast  response  time.  However,  a  severe 
drift  problem  was  encountered  in  using  the  transducer  which 
appeared  insoluble  without  going  into  an  extensive  elec¬ 
tronic  circuitry  development  program.  Consequently,  no 
reproducible  pressure  measurements  were  obtained  with  the 
electrokinetic  transducer.  Because  of  its  high  sensitivity 
and  fast  response,  however,  the  electrokinetic  transducer 
shows  some  promise  as  a  pressure  sensor  if  the  drift  problem 
can  be  solved. 

The  use  of  a  gas  laser  for  measuring  low  velocities  in 
liquids  has  received  considerable  attention  in  recent  years. 
The  laser  velocimeter  has  virtually  an  instantaneous  response 
to  changes  in  flow  velocity.  Limitations  of  time  and 
financial  rescouces  did  not  permit  the  exploration  of  the 
use  of  a  gas  laser  in  the  present  investigation.  Because 
of  the  serious  difficulties  encountered  with  the  micromano¬ 
meter,  however,  it  was  concluded  that  the  use  of  a  micro- 
manometer  for  measuring  small  velocities  in  liquids  is  not 
practical.  The  use  of  the  gas  laser  should  be  considered 
in  future  experiments  of  the  type  described  in  this  inves¬ 
tigation. 

4.3  EXPERIMENTAL  RESULTS 

During  the  course  of  the  experimental  investigation, 
it  became  apparent  that  velocity  profiles  could  not  be 
measured  using  the  micromanometer  because  of  the  slow 
response  and  leak  problem  previously  described.  Consequently, 
the  experimental  effort  was  concentrated  toward  static 
pressure  distribution  measurements  in  the  entrance  region. 

A  limited  number  of  data  points  were  obtained  with  no  mag¬ 
netic  field  using  water  as  the  working  fluid.  The  experi¬ 
mental  results  are  shown  in  Fig.  35. 

A  theoretical  curve  of  entrance  region  pressure  drop 
was  calculated  as  a  function  of  liquid  flow  rate  from  Eqs. 

(60)  through  (62)  and  is  shown  as  the  solid  curve  in  Fig. 

35.  The  agreement  between  the  experimentally  measured 
entrance  region  pressure  drop  and  the  theoretical  curve 
is  considered  as  satisfactory  verification  of  the  theoretical 
model  employed  in  the  present  investigation. 
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V 

SUMMARY 


In  this  investigation  an  analysis  of  the  influence  of 
Hall  current  on  flow  development  in  the  entrance  region  of 
a  parallel  plate  MHD  channel  has  been  presented.  Numerical 
results  have  been  presented  showing  the  influence  of  uyr 
on  the  entrance  length,  the  entrance  region  boundary  layer 
development,  the  entrance  region  pressure  distribution, 
and  the  entrance  region  wall  shear  stress  distribution. 
Constant  electrical  conductivity  was  assumed  in  this  portion 
of  the  analysis.  A  new  linearization  scheme  was  developed 
for  analyzing  MHD  entrance  region  flows  which  has  features 
superior  to  a  conventional  boundary  layer  -  core  flow  model 
of  the  entrance  region.  As  an  alternate  method  of  analyzing 
the  MHD  entrance  region  flow,  the  von-Karman-Pohlhausen 
integral  method  was  generalized  to  accommodate  an  assumed 
polynomial  velocity  profile  of  arbitrary  degree.  Closed 
form  expressions  for  the  velocity  profile  coefficients  were 
obtained,  and  the  range  of  values  of  the  shape  factor  X* 
for  which  the  integral  method  is  valid  was  established  as 
a  function  of  the  degree  of  the  polynomial  assumed. 

In  addition  to  the  entrance  region  results,  two  methods 
of  analyzing  fully  developed  flows  with  arbitrary  variation 
of  conductivity  were  developed.  These  methods  are  the  B.  G. 
Galerkin  method  and  the  Fredholm  integral  equation  method. 
Numerical  results  of  calculated  velocity  profiles  were 
presented  for  several  conductivity  profiles  across  the 
channel.  Although  limitations  of  time  and  financial  re¬ 
sources  did  not  allow  these  methods  to  be  extended  to  the 
entrance  region,  it  is  felt  that  such  an  extension  is 
possible . 

A  new  MHD  channel  flow  facility  was  designed  and  built 
for  obtaining  pressure  distribution  and  velocity  profile 
data  in  the  entrance  region  of  a  parallel  plate  channel. 

The  facility  is  designed  to  use  potassium  chloride  (KC1) 
as  the  working  fluid  with  a  maximum  magnetic  field  strength 
achievable  of  22,000  gauss.  Because  of  serious  instrumen¬ 
tation  problems  encountered  in  measuring  total  velocity 
pressures,  velocity  profile  data  was  not  obtained.  Only 
entrance  region  pressure  distribution  data  was  obtained. 

This  channel  appears  to  be  unique  in  that  no  other 
experimental  facilities  have  been  described  in  the  open 
literature  for  obtaining  MHD  entrance  region  velocity  pro¬ 
files.  The  use  of  a  gas  laser  velocimeter  appears  to  be 
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promising,  and  a  future  experimental  effort  using  a  laser 
for  measuring  the  velocity  profiles  would  represent  a 
significant  contribution  to  a  detailed  understanding  of 
MHD  entrance  region  flows. 
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Fig.  la  MHD  Flow  Geometry 
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Fig.  lb  Geometry  for  Generator  and  Accelerator  Modes 


51 


AEDC-TR-67-79 


A  E  DC-T  R-6?  *79 


A  E  DC-T  R-67-79 


0009 


0008 

0.007 


0.006 

0.005 

0.004 

0.003 


0.002 


2 


2 


AEDC-TR-67-79 


C-TR-67-79 


A  E  DC-T  R  -67-79 


+  ENTRANCE  LENGTH 
-  ASYMPTOTIC  VALUE  OF  K 


Fig.  13  Entrance  Region  Pressure  Correction  for  Ha  =  2 
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Fig.  14  Entrance  Region  Pressure  Correction  for  Ha  =  6 


Fig.  16  Asymptotic  Values  of  K 
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Fig.  18  Entrance  Region  Wall  Shear  Stress  Distribution  for  Ha  =6 
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Fig.  21  Longitudinal  Variation  of  Entrance  Region  Centerline  Velocity 
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Fig,  22  Velocity  Profiles  for  Case  I  with  M  =  3f  'Dy*  =  0 
Using  Fredholm  Integral  Equation  Formulation 


72 


A  ED  C-T  R-67-79 


Fig.  23  Velocity  Profiles  for  Case  II  with  M  =  3,  <ty*  =  0 
Using  Fredholm  Integral  Equation  Formulation 
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Fig.  24  Velocity  Profiles  for  Cose  II  with  M  ~  3,  Ho*  =  2 
Using  Fredholm  Integral  Equation  Formulation 
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Fig.  25  Velocity  Profiles  for  Case  II  with  M  =  3,  Ha*  =  4 
Using  Fredholm  Integral  Equation  Formulation 
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Fig.  26  Velocity  Profile*  for  Cose  II  with  M  =  3,  Ho*  =  8 
Using  Fredholm  Integral  Equation  Formulation 


6 


2 


A  E  DC-T  R  -67*79 


Fig.  28  Velocity  Profiles  for  Case  III  with  M  =  3,  Ha*  =  4 
Using  Fredholm  Integral  Equation  Formulation 
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Fig.  30  Short  Circuit  Velocity  Profiles  for  Cose  IV  by  the  B.G.  Golerkin  Method 
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Fig,  31  Open  Circuit  Velocity  Profiles  for  Case  IV  by  the  B,  G.  GaJerkin  Method 
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Fig.  35  Experimental  Entrance  Region  Pressure  Drop  for  Ha  =0 
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TABLE  2 

EFFECT  OF  THE  DEGREE  OF  THE  POLYNOMIAL  VELOCITY 
PROFILE  ON  THE  ENTRANCE  LENGTH 
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VALUES  OF  OPEN  CIRCUIT  <£y  FOR  CASE  IV 
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